########################################################################
# PLS 206 — Week 5: selection on pure noise
# Grey Monroe, Fall 2026
########################################################################

library(MASS)   # stepAIC(); comes with R

########################################################################
# 1. One run: y and 20 predictors, all random numbers
########################################################################

noise <- data.frame(y = rnorm(100), matrix(rnorm(100 * 20), 100))
head(round(noise[, 1:6], 2))      # nothing is related to anything

step_n <- stepAIC(lm(y ~ ., data = noise), trace = 0)
summary(step_n)                   # any kept predictor with p < 0.05?

# run the three lines above again: a new data set each time

########################################################################
# 2. Repeat 1,000 times: choose on 50 rows, then test on the same 50
#    rows or on the other 50 (about 40 seconds)
########################################################################

one_run <- function() {
  nz <- data.frame(y = rnorm(100), matrix(rnorm(100 * 20), 100))
  a <- nz[1:50, ]                 # rows used to choose
  b <- nz[51:100, ]               # fresh rows, never used to choose
  m <- stepAIC(lm(y ~ ., data = a), trace = 0)
  kept <- attr(terms(m), "term.labels")
  if (length(kept) == 0) return(c(same = FALSE, fresh = FALSE))
  c(same  = any(summary(m)$coefficients[-1, 4] < 0.05),
    fresh = any(summary(lm(reformulate(kept, "y"), data = b))$coefficients[-1, 4] < 0.05))
}

set.seed(1)
runs <- replicate(1000, one_run())
rowMeans(runs)    # share of runs with at least one 'significant' predictor